Preface

This repository contains teaching and learning materials prepared and used during “Introduction to biostatistics and machine learning” course, organized by NBIS, National Bioinformatics Infrastructure Sweden. The course is open for PhD students, postdoctoral researcher and other employees within Swedish universities. The materials are geared towards life scientists wanting to be able to understand and use basic statistical and machine learning methods. More about the course https://uppsala.instructure.com/courses/51998/

Descriptive statistics

Learning outcomes:

• be aware of data types
• compute measures of location, including mean and median
• compute measures of spread, including quantiles, variance and standard deviation
• understand the concept random sample
• compute population mean and variance
• compute sample mean and variance

Data types

Data can be divided into different types; categorical and quantitative (numeric). How to summarize and analyze your data depends on the type.

Categorical data types are divided into;

• Nominal: named. Ex: dead/alive, healthy/sick, WT/mutant, AA/Aa/aa, male/female, red/green/blue
• Ordinal: named and ordered. Ex: pain (weak, moderate, severe), AA/Aa/aa, very young/young/middle age/old/very old, grade I, II, III, IV

Quantitative data types are divided into;

• Discrete: finite or countable infinite values. Ex. counts, number of cells, number of reads
• Continuous: infinitely many uncountable values. Ex. height, weight, concentration

Categorical data

Data summarization and plots

Categorical data can be summarized by counting the number of observations of each category and summarizing in a frequency table or bar plot. Alternatively, the proportions (or percentages) of each category can be calculated.

If you want to follow this example, you can download the data here; mice.csv. You get the subset used in this example by the following commands;

## first read the full data set into R
mice <- read.csv("mice.csv")

## Then extract the specific subset used in this example
m10 <- subset(mice, subset=week==5 & id %in% 1:10, select = c(id, gender, weight))

In this example we have only ten observations (mice) and the full data can actually be shown in a table.

We are interested in the gender distribution in our group of mice. Count the frequency of male/female mice and summarize in a table. Also, the fraction or percentage can be useful.

The frequencies can also be shown in a barplot.

ggplot(m10, aes(x=gender)) + geom_bar()
barplot(table(m10$gender))

(#fig:m10barplot)The number of male and female mice shown in barplots generated using ggplot and basic R graphics.

patients: {L, L, L, R, L, R, R, L, L, R, R, R, R, R, L, R, R, R, R, R, R, R, R, R, R, L, L, R, R, R}

controls: {R, L, R, R, L, R, R, R, R, R, R, R, R, R, R, R, R, R, R, R, R, R, L, R, R, L, R, R, R, R, R, R, R, R, R, R, R, R, R, R}

Summarized as

or in a contingency table;

Data can be summarized in barplots in several ways;

## Using ggplot to create barplots
ggplot(hand, aes(x=group, fill=handedness)) + geom_bar()
ggplot(hand, aes(x=group, fill=handedness)) + geom_bar(position="dodge")
ggplot(hand, aes(x=group, fill=handedness)) + geom_bar(position="fill") + ylab("Fraction")

## Using basic R graphics to create barplots
tab <- table(hand$handedness, hand$group)
barplot(tab)
barplot(tab, beside=TRUE)
tabperc <- tab
tabperc[,1] <- 100*tab[,1]/sum(tab[,1])
tabperc[,2] <- 100*tab[,2]/sum(tab[,2])
barplot(tabperc)

(#fig:barlefthand)Left-handedness in patient and control groups.

Quantitative data

Quantitative data (both discrete and continuous) can be visualized in a histogram;

(#fig:sum10dice)Throw 10 dice and count the total number of dots. Repeat the experiment 1000 times. This histogram summarize the results, i.e. the total number of dots when throwing 10 dice.

or a density plot;

(#fig:dicedens)Density plot over the total number of dots when throwing 10 dice.

(#fig:hist)Histogram over weight of 2000 5 weeks old mice, colored according to gender.

Summary statistics for numeric data are usually divided into measures of location and spread.

Measures of location

• Mode: the most common value, can be computed also for categorical data
• Median: The value that divide the ordered data values into two equally sized groups. 50% of the values are below the median.

• Mean: the arithmetic mean, also called the average

For n onservations x₁, x₂, …, x_n, the mean value is calculated as;

$$\bar x = \frac{x_1+x_2+\dots+x_n}{n} = \frac{1}{n}\displaystyle\sum_{i=1}^n x_i$$ Note, several very different distributions can still have the same mean value.

(#fig:mean35)All these distributions have the same mean value, 3.50.

Measures of spread

• Quartiles - the three values that divide the data values into four equally sized groups.

  • Q1. First quartile. 25% of the values are below Q1. Divides the values below the median into equally sized groups.
  • Q2.´ Second quartile. 50% of the values are below Q2. Q2 is the median.
  • Q3. Third quartile. 75% of the values are below Q3. Divides the values above the median into equally sized groups.

• IQR: interquartile range: Q3 - Q1
• Variance, σ². The variance is the mean squared distance from the mean value.
• Standard deviation, $\sigma = \sqrt{\sigma^2}$.

Variance and standard deviation

The variance of a set of observations is their mean squared distance from the mean value;

$$\sigma^2 = \frac{1}{n} \sum_{i=1}^n (x_i - \bar x)^2.$$ the variance is measured in the square of the unit in which x was measured. a commonly used measured on the same unit as x is the standard deviation, defined as the square root of the variance;

$$\sigma = \sqrt{\frac{1}{n} \sum_{i=1}^n (x_i - \bar x)^2}$$ The denominator n is commonly replaced by n − 1 and the sample standard deviation is calculated instead;

$$s = \sqrt{\frac{1}{n-1} \sum_{i=1}^n (x_i - \bar x)^2}.$$ The latter formula is used if we regard the collection of observations x₁, …, x_n as a sample drawn from a large population of possible observations.

I we want to describe the variance/standard deviation only of our set of observations, the former formula should be used, calculation of a population standard deviation σ (i.e. we consider the set of observations to be the full population).

If instead, we want to estimate the variance of a larger population from which our smaller sample is drawn, we should calculate the sample standard deviation, s.

Exercises: Descriptive statistics

Solutions: Descriptive statistics

@ref(exr:baby) Smokers: 5 (50%) yes baby weight (kg) mean (sd): 3.3 (0.44) gender: 6 (60%) F mother weight (kg) mean(sd): 64 (8.5) mother age mean(sd): 28.5 (5.8) partity mean(sd): 1.5 (1.4) could also be handled as categorical (ordinal) and report frequencies and percentages. married: 4 (40%) yes

Did you caompute standard deviations that were aslightly different? Then you probably computed the sample standard deviation, which could actually be what you want to report. When do you want to compute sample standard deviation?

@ref(exr:cholesterol)

1. $$\bar X = \frac{1}{4}\sum_{i=1}^4 X_i \\ X_i \sim N(\mu, \sigma) \\ \bar X \sim N\left(\mu, \frac{\sigma}{\sqrt{n}}\right) = N(202, 20) $$
2. 0.0019

@ref(exr:pill2)

1. 12.62
2. 0.067
3. 0.26
4. 0.22
5. 0.0028

Probability theory

Learning outcomes

• understand the concept of random variables
• understand the concept of probability
• understand and learn to use resampling to compute probabilities
• understand the concept probability mass function
• understand the concept probability density function
• understand the concept cumulative distribution functions
• use normal distribution
• understand the central limit theorem

Introduction to probability

Some things are more likely to occur than others. Compare:

• the chance of the sun rising tomorrow with the chance that no-one is infected with COVID-19 tomorrow
• the chance of a cold dark winter in Stockholm with the chance of no rainy days over the summer months in Stockholm

We intuitively believe that the chance of sun rising or dark winter occurring are enormously higher than COVID-19 disappearing over night or having no rain over the entire summer. Probability gives us a scale for measuring the likeliness of events to occur. Probability rules enable us to reason about uncertain events. The probability rules are expressed in terms of sets, a well-defined collection of distinct objects.

Suppose we perform an experiment that we do not know the outcome of, i.e. we are uncertain about the outcome. We can however list all the outcomes that might occur.

• sample space is the set S of these possible outcomes of the experiment, e.g. getting 1, 2 etc. on the 6-sided dice S = {1, 2, 3, 4, 5, 6}
• an event is a subset of the sample space
• an event is said to occur if the outcome of the experiment belong to this set
• The complement, E′, of the event E contains all the outcomes in S which are not in E
• Two sets E and F, such that E ∩ F = ∅, are said to be disjoint

Axioms of probability

1. 0 ≤ P(E) ≤ 1 for any event E ⊆ S
2. P(S) = 1
3. if E, F are disjoint events, then P(E∪F) = P(E) + P(F)

Common rules of probability

Based on the axioms the following rules of probability can be proved.

• Complement rule: let E ⊆ S be any event, then P(E′) = 1 − P(E)
• Impossible event: P(∅) = 0
• Probability of a subset: let E, F ⊆ S be events such that E ⊆ F then P(F) ≥ P(E)
• Addition rule: let E, F ⊆ S be any two events, then P(E∪F) = P(E) + P(F) − P(E∩F)

Conditional probability

Let E, F ⊆ S be two events that P(E) > 0 then the conditional probability of F given that E occurs is defined to be: $$P(F|E) = \frac{P(E\cap F)}{P(E)}$$

Product rule follows conditional probability: let E, F ⊆ S be events such that P(E) > 0 then: P(E∩F) = P(F|E)P(E) = P(E|F)P(F)

The urn model

The urn model is a simple model commonly used in statistics and probability. In the urn model, real objects (such as people, mice, cells, genes, molecules, etc) are represented by balls of different colors. A fair coin can be represented by an urn with two balls representing the coins two sides. A group of people can be modelled in an urn model, if age is the variable of interest, we write the age of each person on the balls. If we instead are interested in if the people are allergic to pollen or not, we color the balls according to allergy status.

(#fig:urns)Urn models of a fair coin, age of a group of people, pollen allergy status of a group of people.

In the urn model every unit (ball) is equally likely of being selected. This means that the urn model us well suited to represent flipping a fair coin. However, a biased coin can also be modelled using an urn model, by changing the number of balls that represent each side of the coin.

By drawing balls from the urn with (or without) replacement, probabilities and other properties of the model can be inferred.

Random variables

The outcome of a random experiment can be described by a random variable. Whenever chance is involved in the outcome of an experiment the outcome is a random variable.

A random variable can not be predicted exactly, but the probability of all possible outcomes can be described. The sample space is the set of all possible outcomes of a random variable. Note, the sample space is not always countable.

A random variable is usually denoted by a capital letter, X, Y, Z, …. Values collected in an experiment are observations of the random variable, usually denoted by lowercase letters x, y, z, ….

The population is the collection of all units studied.

A sample is a subset of the population.

Example random variables and probabilites:

• The weight of a random newborn baby, W. P(W>4.0kg)
• The smoking status of a random mother, S. P(S=1)
• The hemoglobin concentration in blood, Hb. P(Hb<125g/L)
• The number of mutations in a gene
• BMI of a random man
• Weight status of a random man (underweight, normal weight, overweight, obese)
• The result of throwing a die

Conditional probability can be written for example P(W≥3.5|S=1), which is the probability that X ≥ 3.5 if S = 1, in words “the probability that a smoking mother has a baby with birth weight of 3.5 kg or more”.

Discrete random variables

A discrete random number has countable number of outcome values, such as {1,2,3,4,5,6}; {red, blue, green}; {tiny, small, average, large, huge} or all integers.

A discrete random variable can be described by its probability mass function, pmf.

The probability that the random variable, X, takes the value x is denoted P(X=x) = p(x). Note that:

1. 0 ≤ p(x) ≤ 1, a probability is always between 0 and 1.
2. ∑p(x) = 1, the sum over all possible outcomes is 1.

The probability mass function of a fair six-sided die can be summarized in a table;

or in a barplot;

(#fig:die)Probability mass function of a die.

(#fig:CFU)Probability mass distribution of the number of bacterial colonies on an agar plate.

Expected value

When the probability mass function is know the expected value of the random variable can be computed.

$$E[X] = \mu = \sum_{i=1}^N x_i p(x_i)$$ For a uniform distribution, where every object has the same probability (in the urn model, every object is represented by one ball), the expected value can be computed as the sum of all objects divided by the total number of objects;

$$E[X] = \mu = \frac{1}{N}\sum_{i=1}^N x_i$$ Linear transformations and combinations

E(aX) = aE(X)

E(X+Y) = E(X) + E(Y)

E[aX+bY] = aE[X] + bE[Y]

Variance

The variance is a measure of spread and is defined as the expected value of the squared distance from the expected value;

$$var(X) = \sigma^2 = E[(X-\mu)^2] = \sum_{i=1}^n (x_i-\mu)^2 p(x_i)$$ Linear transformations and combinations

var(aX) = a²var(X)

For independent random variables X and Y

var(aX+bY) = a²var(X) + b²var(Y)

Simulate distributions

Once the distribution is known, we can compute probabilities, such as P(X=x), P(X<x) and P(X≥x). If the distribution is not known, simulation might be the solution.

The outcome of a single coin toss is a random variable, X with two possible outcomes {H, T}. We know that P(X=H) = 0.5. The random variable of interest is the number of heads in 20 coin tosses, Y. The probability that we need to compute is P(Y≥15).

(#fig:coinurn)A coin toss. Urn model with one black ball (heads) and one white ball (tails).

A single coin toss can be modelled by an urn with two balls. When a ball is drawn randomly from the urn, the probability to get the black ball (heads) is P(X=H) = 0.5.

If we want to simulate tossing 20 coins (or one coin 20 times) we can use the same urn model, if the ball is replaced after each draw.

In R we can simulate random draws from an urn model using the function sample.

# A single coin toss
sample(c("H", "T"), size=1)

## [1] "H"

# Another coin toss
sample(c("H", "T"), size=1)

## [1] "H"

Every time you run sample a new coin toss is simulated.

The argument size tells the function how many balls we want to draw from the urn. To draw 20 balls from the urn, set size=20, remember to replace the ball after each draw!

# 20 independent coin tosses
(coins <- sample(c("H", "T"), size=20, replace=TRUE))

##  [1] "H" "H" "H" "H" "T" "H" "H" "T" "T" "T" "H" "T" "H" "H" "T" "H" "T" "T" "T"
## [20] "T"

How many heads did we get in the 20 random draws?

# How many heads?
sum(coins == "H")

## [1] 10

We can repeat this experiment (toss 20 coins and count the number of heads) several times to estimate the distribution of number of heads in 20 coin tosses.

To do the same thing several times we use the function replicate.

To simulate tossing 20 coins and counting the number of heads 10000 times, do the following;

Nheads <- replicate(10000, {
  coins <- sample(c("H", "T"), size=20, replace=TRUE)
  sum(coins == "H")
})

Plot distribution of the number of heads in a histogram.

hist(Nheads, breaks=0:20-0.5)

Now, let’s get back to the question; when tossing 20 coins, what is the probability of at least 15 heads?

P(X≥15)

Count how many times out of our 10000 exeriments the number is 15 or greater

sum(Nheads >= 15)

## [1] 200

From this we conclude that

P(X≥15)= 200/10000 = 0.02

Parametric discrete distributions

Bernoulli trial

A Bernoulli trial is a random experiment with two outcomes; success and failure. The probability of success, P(success) = p, is constant. The probability of failure is P(failure) = 1 − p.

When coding it is convenient to code success as 1 and failure as 0.

The outcome of a Bernoulli trial is a discrete random variable, X.

$$p(x) = \left\{ \begin{array}{ll} p & \mathrm{if}\,x=1\mathrm,\,success\\ 1-p & \mathrm{if}\,x=0\mathrm,\,failure \end{array} \right.$$

Using the definitions of expected value and variance it can be shown that;

$$E[X] = p\\ var(X) = p(1-p)$$

Binomial distribution

The number of successes in a series of independent and identical Bernoulli trials is a discrete random variable, X.

$X = \sum_{i=0}^n Z_i,$

where all Z_i describe the outcome of independent and identical Bernoilli trials with probability p for success (P(Z_i=1) = p).

The probability mass function of X is called the binomial distribution. In short we use the notation;

X ∈ Bin(n,p)

The probability mass function is

$$P(X=k) = {n \choose k} p^k (1-p)^{n-k}$$ It can be shown that

$$E[X] = np\\ var(X) = np(1-p)$$

The binomial distribution occurs when sampling n objects with replacement from an urn with objects of two types, of which the interesting type has probability p.

The probability mass function, P(X=k) can be computed using the R function dbinom and the cumulative distribution function P(X≤k) can be computed using pbinom.

Hypergeometric distribution

The hypergeometric distribution occurs when sampling n objects without replacement from an urn with N objects of two types, of which the interesting type has probability p.

The probability density function

$$P(X=k) = \frac{{Np \choose x} {N-Np \choose n-x}}{N \choose n}$$ can be computed in R using dhyper and the cumulative distribution function P(X≤k) can be computed using phyper.

Poisson distribution

The Poisson distribution describe the number of times a rare event occurs in a large number of trials.

A rare disease has a very low probability for a single individual. The number of individuals in a large population that catch the disease in a certain time period can be modelled using the Poisson distribution.

The probability mass function;

$$P(X=k) = \frac{\mu}{k!}e^{-\mu},$$ where μ is the expected value, which is μ = nπ, where n is the number of objects sampled from the population and π is the probability of a single object.

The Poisson distribution can approximate the binomial distribution if n is large (n > 10) and π is small (π < 0.1).

Distributions in R

Probability mass functions, P(X=x), for the binomial, hypergeometric and Poisson distributions can in R can be computed using functions dbinom, dhyper, and dpois, respectively.

Cumulative distribution functions, P(X≤x) can be computed using pbinom, phyper and ppois.

Also, functions for computing an x such that P(X≤x) = q, where q is a probability of interest are available using qbinom, qhyper, and qpois.

Exercises: Discrete random variables

Introduction to probability

Simulation

Parametric discrete distributions

Conditional probability

1. the probability of a positive test result from a person with cancer?

2. the probability of a negative test result from a person without cancer?

3. the probability of having cancer, if the test is positive?

4. the probability of not having cancer, if the test is negative?

5. Connect the four computed probabilities with the following four terms;

• Sensitivity
• Specificity
• Positive predictive value (PPV)
• Negative predictive value (NPV)

Solutions: Discrete random variables

@ref(exr:probcoin)

1. 0.5
2. 0.5

@ref(exr:probdie)

1. $\frac{1}{6}$
2. $\frac{3}{6} = \frac{1}{2}$
3. $\frac{4}{6} = \frac{2}{3}$
4. $1*\frac{1}{6} + 2*\frac{1}{6} + 3*\frac{1}{6} + 4*\frac{1}{6} + 5*\frac{1}{6} + 6*\frac{1}{6} = 3.5$

@ref(exr:Cointoss)

Simulate as in the lecture;

# A single coin toss
sample(c("H", "T"), size=1)

## [1] "T"

# Another coin toss
sample(c("H", "T"), size=1)

## [1] "H"

# 20 independent coin tosses
(coins <- sample(c("H", "T"), size=20, replace=TRUE))

##  [1] "H" "H" "T" "T" "T" "T" "H" "T" "H" "H" "H" "T" "T" "H" "T" "T" "H" "H" "T"
## [20] "H"

# How many heads did we get in these particular 20 draws?
sum(coins == "H")

## [1] 10

## The simulation is about repeating this (20 random draws and summing up the number of heads) many times. To do it 10000 times;
Nheads <- replicate(10000, {
  coins <- sample(c("H", "T"), size=20, replace=TRUE)
  sum(coins == "H")
})

1. Probability of exactly 15 heads

## Numer of times of the 10000 with exactly 15 heads
sum(Nheads==15)

## [1] 156

## divide by 10000 to get the probability
sum(Nheads==15)/10000

## [1] 0.016

## or compute using mean (why does this work?)
mean(Nheads==15)

## [1] 0.016

2. Probability of less than 7 heads

mean(Nheads<7)

## [1] 0.055

3. What is the most probable number of heads?

## plot the distribution and read the graph
hist(Nheads, breaks=0:20-0.5)

## or tabulate
table(Nheads)

## Nheads
##    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18 
##    8   37  145  359  718 1188 1662 1751 1616 1211  710  387  156   41    9    2

4. What is the probability of 5 tails or less?

To get five or less tails out of 20 throws is equal to getting 15 or more heads out of 20.

## probability of 15 heads or more
mean(Nheads>=15)

## [1] 0.021

5. what is the probability of 2 heads or less?

mean(Nheads<=2)

## [1] 0

sum(Nheads<=2)

## [1] 0

## with this low number of observations, more repeats is required to get a more accurate answer

Nheads <- replicate(1000000, {
  coins <- sample(c("H", "T"), size=20, replace=TRUE)
  sum(coins == "H")
})
sum(Nheads<=2)

## [1] 215

mean(Nheads<=2)

## [1] 0.00021

1. 0.015
2. 0.058
3. 10
4. 0.020
5. 0.0002

@ref(exr:Dice)

1. X, the number of sixes when rolling 10 dice.
2. 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10

N <- replicate(100000, sum(sample(1:6, size=10, replace=TRUE)==6))
table(N)

## N
##     0     1     2     3     4     5     6     7     8 
## 16141 32153 29094 15502  5478  1366   234    30     2

hist(N, breaks=(0:11)-0.5)

6. What is the probability to get exactly 2 sixes?
7. On average how many sixes do you get when throwing ten dice?
8. 0.015

## [1] 1632

## [1] 0.016

## [1] 0.016

5. 1
6. 0.29

mean(N==2)

## [1] 0.29

7. 1.7

mean(N)

## [1] 1.7

10*1/6

## [1] 1.7

@ref(exr:Cards)

0.00049

N <- replicate(100000, sum(sample(rep(0:1, c(39,13)), size=5)))
hist(N, breaks=(0:6)-.5)

# P(N==5)
mean(N==5)

## [1] 0.00039

@ref(exr:Pollen)

## Solution using 100 replicates
x <- replicate(100, sum(sample(c(0,0,0,0,0,0,0,1,1,1), size=3, replace=TRUE)))
table(x)

## x
##  0  1  2  3 
## 42 43 13  2

mean(x==0)

## [1] 0.42

## Solution using 1000 replicates
x <- replicate(1000, sum(sample(c(0,0,0,0,0,0,0,1,1,1), size=3, replace=TRUE)))
table(x)

## x
##   0   1   2   3 
## 358 429 193  20

mean(x==0)

## [1] 0.36

## Solution using 100000 replicates
x <- replicate(100000, sum(sample(c(0,0,0,0,0,0,0,1,1,1), size=3, replace=TRUE)))
table(x)

## x
##     0     1     2     3 
## 34253 44113 18927  2707

mean(x==0)

## [1] 0.34

@ref(exr:Pollen2)

## Solution using 100000 replicates
x <- replicate(100000, sum(sample(rep(c(0, 1), c(14, 6)), size=3, replace=FALSE)))
table(x)

## x
##     0     1     2     3 
## 32080 47809 18368  1743

mean(x==0)

## [1] 0.32

@ref(exr:Pollen3)

## Solution using 100000 replicates
x <- replicate(100000, sum(sample(rep(c(0, 1), c(140, 60)), size=3, replace=FALSE)))
table(x)

## x
##     0     1     2     3 
## 34182 44396 18785  2637

mean(x==0)

## [1] 0.34

Parametric discrete distributions

@ref(exr:pollenparam)

## 1.6 Solution using the Binomial distribution
pbinom(0, 3, 0.3)

## [1] 0.34

## 1.7 Solution using the hypergeometric distribution
phyper(0, 6, 20-6, 3)

## [1] 0.32

## 1.8 Solution using the hypergeometric distribution
phyper(0, 66, 200-60, 3)

## [1] 0.31

@ref(exr:GSEA)

phyper(20, 200, 20000-200, 1000, lower.tail=FALSE)

## [1] 0.0011

Conditional probability

@ref(exr:diagnostictests)

1. P(pos|cancer) = 16/20, sensitivity
2. P(neg|notcancer) = 882/980, specificity
3. P(cancer|pos) = 16/114, PPV
4. P(notcancer|neg) = 882/886, NPV

Continuous random variable

A continuous random number is not limited to discrete values, but any continuous number within one or several ranges is possible.

Examples: weight, height, speed, intensity, …

A continuous random variable can be described by its probability density function, pdf.

(#fig:unnamed-chunk-32)Probability density function of the weight of a newborn baby.

The probability density function, f(x), is defined such that the total area under the curve is 1.

∫_−∞^∞f(x)dx = 1

The area under the curve from a to b is the probability that the random variable X takes a value between a and b.

P(a≤X≤b) = ∫_a^bf(x)dx

The cumulative distribution function, cdf, sometimes called just the distribution function, F(x), is defined as:

F(x) = P(X≤x) = ∫_−∞^xf(x)dx

P(X≤x) = F(x)

As we know that the total probability (over all x) is 1, we can conclude that

P(X>x) = 1 − F(x) and thus

P(a<X≤b) = F(b) − F(a)

Parametric continuous distributions

Two important parameters of a distribution is the expected value, μ, that describe the distributions location and the variance, σ², that describe the spread.

The expected value, or population mean, is defined as;

E[X] = μ = ∫_−∞^∞xf(x)dx We will learn more about the expected value and how to estimate a population mean from a sample later in the course.

The variance is defined as the expected value of the squared distance from the population mean;

σ² = E[(X−μ)²] = ∫_−∞^∞(x−μ)²f(x)dx

The square root of the variance is called the standard deviation, σ.

Normal distribution

The normal distribution (sometimes referred to as the Gaussian distribution) is a common probability distribution and many continuous random variables can be described by the normal distribution or be approximated by the normal distribution.

The normal probability density function

$$f(x) = \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{1}{2} \left(\frac{x-\mu}{\sigma}\right)^2}$$

describes the distribution of a normal random variable, X, with expected value μ and standard deviation σ, e and π are two common mathematical constants, e ≈ 2.71828 and π ≈ 3.14159.

In short we write X ∼ N(μ,σ).

The bell-shaped normal distributions is symmetric around μ and f(x) → 0 as x → ∞ and as x →  − ∞.

As f(x) is well defined, values for the cumulative distribution function F(x) = ∫_−∞^xf(x)dx can be computed.

If X is normally distributed with expected value μ and standard deviation σ we write:

X ∼ N(μ,σ)

Using transformation rules we can define

$$Z = \frac{X-\mu}{\sigma}, \, Z \sim N(0,1)$$

Values for the cumulative standard normal distribution, F(z), are tabulated and easy to compute in R using the function pnorm.

(#fig:FZ)The shaded area under hte curve is the tabulated value P(Z≤z) = F(z).

Some value of particular interest:

$$F(1.64) = 0.95\\ F(1.96) = 0.975$$

As the normal distribution is symmetric F(-z) = 1 - F(z)

$$F(-1.64) = 0.05\\ F(-1.96) = 0.025$$

P(−1.96<Z<1.96) = 0.95

Sum of two normal random variables

If X ∼ N(μ₁,σ₁) and Y ∼ N(μ₂,σ₂) are two independent normal random variables, then their sum is also a random variable:

$$X + Y \sim N(\mu_1 + \mu_2, \sqrt{\sigma_1^2 + \sigma_2^2})$$

and

$$X - Y \sim N(\mu_1 - \mu_2, \sqrt{\sigma_1^2 + \sigma_2^2})$$ This can be extended to the case with n independent and identically distributed random varibles X_i (i = 1…n). If all X_i are normally distributed with mean μ and standard deviation σ, X_i ∈ N(μ,σ), then the sum of all n random variables will also be normally distributed with mean nμ and standard deviation $\sqrt{n} \sigma$.

Central limit theorem

As a result of central limit theorem, the distribution of fractions or mean values of a sample follow the normal distribution, at least if the sample is large enough (a rule of thumb is that the sample size n > 30).

##Population mean
(mu <- mean(fat$BMI))

## [1] 25

##Population variance
(sigma2 <- var(fat$BMI)/nrow(fat)*(nrow(fat)-1))

## [1] 13

##Population standard variance
(sigma <- sqrt(sigma2))

## [1] 3.6

Randomly sample 3, 5, 10, 15, 20, 30 men and compute the mean value, m. Repeat many times to get the distribution of mean values.

Note, mean is just the sum divided by the number of samples n.

χ²-distribution

The random variable $Y = \sum_{i=1}^n X_i^2$ is χ² distributed with n − 1 degrees of freedom, if X_i are independent identically distributed random variables X_i ∈ N(0,1).

In short Y ∈ χ²(n−1).

(#fig:Xdistr)The χ²-distribution.

F-distribution

The ratio of two χ²-distributed variables divided by their degrees of freedom is F-distributed

(#fig:Fdistr)The F-distribution

t-distribution

The ratio of a normally distributed variable and a χ²-distributed variable is t-distributed.

(#fig:exampletdistr)The t-distribution.

Distributions in R

Probability density functions for the normal, t, χ² and F distributions can in R can be computed using functions dnorm, dt, dchisq, and df, respectively.

Cumulative distribution functions can be computed using pnorm, pt, pchisq and pf.

Also, functions for computing an x such that P(X<x) = q, where q is a probability of interest are available using qnorm, qt, qchisq and qf.

Exercises: Continuous random variables

Solutions

@ref(exr:normtable)

1. P(Z<1.64) = P(Z≤1.64) = [fromtable] = 0.9495
2. P(Z>−1.64) = [symmetry] = P(Z<1.64) = 0.9495
3. P(−1.96<Z) = [symmetry] = P(Z<1.96) = [fromtable] = 0.975
4. P(Z<2.36) = [fromtable] = 0.9909
5. Standard normal table gives that a = 1.64
6. P(Z>b) = 0.975, symmetry gives that P(Z<−b) = 0.975. Look-up in the standard normal table gives that  − b = 1.96, hence b =  − 1.96.

@ref(exr:ztransform)

1. $P(X<5) = P\left(\frac{X-3}{2}\right < \frac{5-3}{2}) = P(Z<1) = [±textrm{from table}] = 0.8413$
2. $P(3<X<5) = P\left(\frac{3-3}{2} < \frac{X-3}{2} < \frac{5-3}{2}\right) = P(0<Z<1) = P(Z<1) - P(Z<0) = 0.8413 - 0.5000 = 0.3413$
3. $P(X> \geq 7) = P\left(\frac{X-3}{2} \geq \frac{7-3}{2}) = P(Z \geq 2) = 1 - P(Z < 2) = 1 - 0.9772 = 0.0228$

@ref(exr:sumdistr)

1. $P(Hb<158) = P(Z<\frac{158-188}{14}) = P(Z<-2.14) = [table] = 0.016$ or use R

pnorm(158, mean=188, sd=14)

## [1] 0.016

2. Probability of one man not being anemic; 1 − 0.016 = 0.984. Probability of 10 selected men not being anemic (binomial distribution) 0.984¹⁰ = 0.95

@ref(exr:poisson) a)

n <- 10000
p <- 3/100000
ppois(0, n*p)

## [1] 0.74

ppois(1, n*p, lower.tail=FALSE)

## [1] 0.037

Random sample

In many (most) experiments it is not feasible (or even possible) to examine the entire population. Instead we study a random sample.

A random sample is a random subset of individuals from a population.

A simple random sample is a random subset of individuals from a population, where every individual has the same probability of being choosen.

The urn model to perform simple random sampling

Let every individual in the population be represented by a ball. The value on each ball is the measurement we are interested in, for example height, shoe size, hair color, healthy/sick, type of cancer/no cancer, blood glucose value, etc.

Draw n balls from the urn, without replacement, to get a random sample of size n.

Sample properties

Summary statistics can be computed for a sample, such as the sum, proportion, mean and variance.

Notation: A random sample x₁, x₂, …, x_n from a distribution, D, consists of n observations of the independent random variables X₁, X₂, …, X_n all with the distribution D.

Sample proportion

The proportion of a population with a particular property is π.

The number of individuals with the property in a simple random sample of size n is a random variable X. The proportion of individuals in a sample with the property is also a random variable;

$$P = \frac{X}{n}$$ with expected value $$E[P] = \frac{E[X]}{n} = \frac{n\pi}{n} = \pi$$

Sample mean and standard deviation

For a particular sample of size n; x₁, …, x_n, the sample mean is denoted m = x̄. The sample mean is calculated as;

$$m = \bar x = \frac{1}{n}\displaystyle\sum_{i=1}^n x_i$$ and the sample variance as;

$$s^2 = \frac{1}{n-1} \sum_{i=1}^n (x_i-m)^2$$

Note that the mean of n independent identically distributed random variables, X_i is itself a random variable;

$$\bar X = \frac{1}{n}\sum_{i=1}^n X_i,$$ If X_i ∼ N(μ,σ) then $\bar X \sim N\left(\mu, \frac{\sigma}{\sqrt{n}}\right)$.

When we only have a sample of size n, the sample mean m is our best estimate of the population mean. It is possible to show that the sample mean is an unbiased estimate of the population mean, i.e. the average (over many size n samples) of the sample mean is μ.

$$E[\bar X] = \frac{1}{n} n E[X] = E[X] = \mu$$

Similarly, the sample variance is an unbiased estimate of the population variance.

Standard error

Eventhough the sample can be used to calculate unbiased estimates of the population value, the sample estimate will not be perfect. The standard deviation of the sampling distribution (the distribution of sample estimates) is called the standard error.

For the sample mean, X̄, the variance is

$$E[(\bar X - \mu)] = var(\bar X) = var(\frac{1}{n}\sum_i X_i) = \frac{1}{n^2} \sum_i var(X_i) = \frac{1}{n^2} n var(X) = \frac{\sigma^2}{n}$$ The standard error of the mean is thus;

$$SEM = \frac{\sigma}{\sqrt{n}}$$ Replacing σ with the sample standard deviation, s, we get an estimate of the standard deviation of the mean;

$$SEM \approx \frac{s}{\sqrt{n}}$$ An alternative definition of standard error of the mean is actually

$$SEM = \frac{s}{\sqrt{n}}$$

Statistical Inference

Learning outcomes:

• to define null and alternative hypothesis
• to perform a hypothesis test using resampling
• to perform a t-test
• to understand and define sampling distribution and standard error
• to compute standard error of mean and proportions
• to compute confidence interval of mean and proportions using the normal approximation
• to compute confidence interval of mean using the t-distribution

Statistical inference is to draw conclusions regarding properties of a population based on observations of a random sample from the population.

Hypothesis test

To perform a hypothesis test is to evaluate a hypothesis based on a random sample.

Typicaly, the hypotheses that are tested are assumptions about properties of the population, such as proportion, mean, mean difference, variance etc.

The null and alternative hypothesis

There are two hypotheses involved in a hypothesis test, the null hypothesis, H₀, and the alternative hypothesis, H₁.

The null hypothesis is in general neutral, “no change”, “no difference between groups”, “no association”. In general we want to show that H₀ is false.

The alternative hypothesis expresses what the researcher is interested in “the treatment has an effect”, “there is a difference between groups”, “there is an association”. The alternative hypothesis can also be directional “the treatment has a positive effect”.

To perform a hypothesis test

1. Define H₀ and H₁
2. Select appropriate test statistic, T, and compute the observed value, t_obs
3. Assume that the H₀ is true and compute the sampling distribution of T.
4. Select an appropriate significance level, α
5. Compare the observed value, t_obs, with the computed sampling distribution under H₀ and compute a p-value. The p-value is the probability of observing a value at least as extreme as the observed value, if H₀ is true.
6. Based on the p-value either accept or reject H₀.

(#fig:examplenull)A null distribution

(#fig:examplepval)The p-value is the probability to observe x_obs or something more extreme, if the null hypothesis is true.

Significance level and error types

The significance level, α = P(false alarm) = P(Reject H₀|H₀ is true).

The significance level is the risk to of false alarm, i.e. to say “I have a hit”, “I found a difference”, when the the null hypothesis (“there is no difference”) is true. The risk of false alarm is control by setting the significance level to a disired value. We do want to keep the risk of false alarm (type I error) low, but at the same time we don’t want to many missed hits (type II error).

The significance level should be set before the hypothesis test is performed. Common values to use are 0.05 or 0.01.

If the p-value is above the significance level, H₀ is accepted.

If the p-value is below the significance level, H₀ is rejected.

Hypothesis test, simulation examples

In these examples the significance level is set to 0.05.

Null and alternative hypotheses

H₀: The proportion of pollen allergy in Uppsala is the same as in Sweden as a whole.

H₁: The proportion of pollen allergy in Uppsala is not the same as in Sweden as a whole.

or expressed differently;

H₀: π = π₀

H₁: π > π₀ where π is the unknown proportion of pollen allergy in the Uppsala population that. π₀ = 0.3 is the proportion of pollen allergy in Sweden.

Test statistic

Here we are interested in the proportion of pollen allergic in Uppsala. An appropriate test statistic could be the number of pollen allergic in a sample of size n = 100, X. As an alternative we can use the proportion of pollen allergic in a sample of size n,

$$P = \frac{X}{n}$$

Let’s use P as our test statistic and compute the observed value, p_obs. In our sample of 100 people from Uppsala the proportion allergic to pollen is p = 42/100 = 0.42.

Null distribution

The sampling distribution of P under H₀ (i.e. when the null hypothesis is true) is what we call the null distribution.

H₀ state that π = 0.3. We can model this using an urn model as follows;

(#fig:pollenurn)An urn model of the null hypothesis π = 0.3. The black balls represent allergic and the white balls non-allergic.

Using this model, we can simulate taking a sample of size 100 many times.

## Urn
rep(c(0, 1), c(7, 3))

##  [1] 0 0 0 0 0 0 0 1 1 1

## Sample 100 times with replacement
sample(rep(c(0, 1), c(7, 3)), 100, replace=TRUE)

##   [1] 0 0 0 1 0 0 1 1 1 0 0 1 0 0 0 0 1 0 1 1 0 1 0 1 0 0 0 0 1 0 1 1 0 0 0 0 1
##  [38] 0 1 0 0 0 0 1 0 0 0 0 1 1 0 0 0 1 0 0 1 0 1 0 1 0 0 1 1 0 0 0 1 1 1 0 0 0
##  [75] 0 0 0 0 0 0 0 1 1 0 0 0 1 1 1 0 1 0 1 1 1 0 1 0 0 0

## Compute proportion of samples that are allergic (1)
sum(sample(rep(c(0, 1), c(7, 3)), 100, replace=TRUE))/100

## [1] 0.3

## Draw samples of size 100 and compute proporion allergic 100000 times
p <- replicate(100000, mean(sample(rep(c(0, 1), c(7, 3)), 100, replace=TRUE)))

Finally plot the distribution

(#fig:pollensampledistr)The sampling distribution.

Compute p-value

Compare the observed value, p_obs = 0.42 to the null distribution.

(#fig:unnamed-chunk-41)The sampling distribution. The observed value is marked by a red vertical line.

The p-value is the probability of getting the observed value or higher, if the null hypothesis is true.

Use the null distribution to calculate the p-value, P(P≥0.42|H₀).

## How many times 
sum(p >= 0.42)

## [1] 689

## p-value
sum(p >= 0.42)/length(p)

## [1] 0.00689

p = P(P≥0.42|H₀) = 0.00689

Accept or reject H₀?

The observed values, mouse weights in grams, are summarized below;

Null and alternative hypotheses

$$ \begin{aligned} H_0: \mu_2 = \mu_1 \iff \mu_2 - \mu_1 = 0\\ H_1: \mu_2>\mu_1 \iff \mu_2-\mu_1 > 0 \end{aligned} $$

where μ₂ is the (unknown) mean body weight of the high-fat mouse population and μ₁ is the mean body-weight of the control mouse population.

Studied population: Female mice that can be ordered from a lab.

Test statistic

Here we are interested in the mean difference between high-fat and control mice.

Mean weight of 12 (randomly selected) mice on ordinary diet, X̄₁. E[X̄₁] = E[X₁] = μ₁

Mean weight of 12 (randomly selected) mice on high-fat diet, X̄₂. E[X̄₂] = E[X₂] = μ₂

The mean difference is also a random variable: D = X̄₂ − X̄₁

Observed values;

## 12 HF mice
xHF <- c(25, 30, 23, 18, 31, 24, 39, 26, 36, 29, 23, 32)
## 12 control mice
xN <- c(27, 25, 22, 23, 25, 37, 24, 26, 21, 26, 30, 24)

##Compute mean body weights of the two samples
mHF <- mean(xHF)
mN <- mean(xN) 
## Compute mean difference
dobs <- mHF - mN

Mean weight of control mice (ordinary diet): x̄₁ = 25.83

Mean weight of mice on high-fat diet: x̄₂ = 28.00

Difference in mean weights: d_obs = x̄₂ − x̄₁ = 2.17

Null distribution

If high-fat diet has no effect, i.e. if H₀ was true, the result would be as if all mice were given the same diet. What can we expect if all mice are fed with the same type of food?

This can be accomplished using permutation

The 24 mice were initially from the same population, depending on how the mice are randomly assigned to high-fat and normal group, the mean weights would differ, even if the two groups were treated the same.

Assume H₀ is true, i.e. assume all mice are equivalent and

1. Randomly reassign 12 of the 24 mice to ‘high-fat’ and the remaining 12 to ‘control’.
2. Compute difference in mean weights

If we repeat 1-2 many times we get the sampling distribution when H₀ is true, the so called null distribution, of difference in mean weights.

## All 24 body weights in a vector
x <- c(xHF, xN)
## Mean difference
(dobs <- mean(x[1:12]) - mean(x[13:24]))

## [1] 2.167

## Permute once
(y <- sample(x))

##  [1] 29 36 39 31 26 24 25 23 25 27 32 37 30 23 30 18 24 23 24 26 22 26 21 25

##Compute mean difference
mean(y[1:12]) - mean(y[13:24])

## [1] 5.167

dnull.perm <- replicate(n = 100000, {
  y <- sample(x)
  ##Mean difference
  mean(y[1:12]) - mean(y[13:24])
})
ggplot(data.frame(d=dnull.perm), aes(x=d)) +
  geom_histogram(bins=25, color="white") +
  theme_bw() +
  geom_vline(xintercept=dobs, color="red")

##Alternatively plot using hist
## hist(dnull.perm)

Compute p-value

What is the probability to get an at least as extreme mean difference as our observed value, d_obs, if H₀ was true?

## Compute the p-value
sum(dnull.perm>dobs)/length(dnull.perm)

## [1] 0.1492

mean(dnull.perm>dobs)

## [1] 0.1492

$P(X_2 - X_1 d_{obs} | H_0) = $ 0.168

Conclusion?

Exercises: Hypothesis tests, resampling

ctrl <- c(197, 186, 157, 170, 193, 188, 175, 186, 177, 191, 168, 193, 191, 189, 188, 192, 179, 186, 197, 203)
iron <- c(187, 218, 196, 210, 206, 178, 181, 193, 172, 202, 169, 221, 183, 222, 185, 174, 192, 192, 162, 211)

Perform a hypothesis test to investigate if the Hb level is affected (increased or decreased) by the iron-rich diet.

Parametric tests

In previous chapters we have computed the sampling distribution using resampling techniques to be able to perform hypothesis tests or compute interval estimates. If the null distribution was already known (or could be computed based on a few assumptions) resampling would not be necessary.

We can follow the same steps as before to perform a hypothesis test:

1. Define H₀ and H₁
2. Select appropriate test statistic, T, and compute the observed value, t_obs
3. Assume that the H₀ and derive the null distribution of the test statistic based on appropriate assumptions.
4. Select an appropriate significance level, α
5. Compare the observed value, t_obs, with the null distribution and compute a p-value. The p-value is the probability of observing a value at least as extreme as the observed value, if H₀ is true.
6. Based on the p-value either accept or reject H₀.

One sample, proportions

The number of allergic individuals in a sample of size n is X and the proportion of allergic persons is P = X/n. X is binomially distributed, but here we can use the Central limit theorem, see @ref(thm:CLT).

As a result of the central limit theorem, the distribution of number or proportion of allergic individuals in a sample of size n is approximately normal. At least if the sample is large enough. A rule of thumb is that the sample size should be n > 30.

Here, the sample size is 100!

The normal distribution has two parameters, mean and standard deviation.

From the binomial distribution we know that E[X] = π and var(X) = nπ(1−π). Hence E[P] = π and $var(P) = \frac{\pi(1-\pi)}{n}$.

The standard error is thus

$$SE=\sqrt{\frac{\pi(1-\pi)}{n}}$$

When the null hypothesis is true π is known and π = 0.3.

Actually these calculations are true in general when a proportion in one sample is compared to a known value.

$$H_0: \pi=\pi_0 \\ H_1: \pi>\pi_0 $$

Ather potential alternative hypothesis are H₁ : π < π₀ or H₁ : π ≠ π₀, but in this particular example we are only interested in the alternative that π > π₀.

If H₀ is true π = π₀ and

$$P \sim N\left(\pi_0, \sqrt{\frac{\pi_0(1-\pi_0)}{n}}\right)$$ An appropriate test statistic is

$$Z = \frac{P-\pi_0}{\sqrt{\frac{\pi_0(1-\pi_0)}{n}}}$$

Z ∈ N(0,1) which makes probabilities easy to compute.

Back to our example, replace P with our observed value p = 0.42 and π₀ = 0.3 and compute our observed

$$Z_{obs} = \frac{0.42-0.3}{\sqrt{\frac{0.3(1-0.3)}{100}}} = 2.62$$

The p-value is the probability of the observed value, or something more extreme, if the null hypothesis is true. If the computed probability is below α = 0.05 our significance threshold, H₀ will be rejected.

p = P(P≥π₀) = P(Z≥Z_obs) = P(Z≥2.62) = 1 − P(Z≤2.62) = [table] = 1 − 0.996 = 0.0044

As 0.0044<0.05 we reject H₀ and conclude that there is reason to believe that the proportion of allergic in Uppsala is greater than 0.3.

One sample, mean

A one sample test of means compares the mean of a sample to a prespecified value.

For example, we might know that the weight of a mouse on normal diet is normally distributed with mean 24.0 g and standard deviation 3 g and want to compare the weight of a sample of 10 mice on high-fat diet to the known mean value for mice on normal diet.

The hypotheses:

$$H_0: \mu = \mu_0 \\ H_1: \mu \neq \mu_0$$

The alternative hypothesis, H₁, above is for the two sided hypothesis test. Other options are the one sided H₁; H₁ : μ > μ₀ or H₁ : μ < μ₀.

If X ∼ N(μ,σ) (this could for example be the weight of a mouse on high-fat diet) then the sample mean $$\bar X \sim N\left(\mu, \frac{\sigma}{\sqrt{n}}\right).$$

If σ is known the test statistic

$$Z = \frac{\bar X - \mu}{\frac{\sigma}{\sqrt{n}}}$$ is normally distributed,  ∼ N(0,1).

For small n and unknown σ, the test statistic

$$t = \frac{\bar X - \mu}{\frac{s}{\sqrt{n}}}$$

is t-distributed with df = n − 1 degrees of freedom.

Two samples, proportions

$$H_0: \pi_1 - \pi_2 = 0\\ H_1: \pi_1 - \pi_2 \neq 0$$

Alternatively, a one sided alternative hypothesis can be used; H₁ : π₁ − π₂ > 0 or H₁ : π₁ − π₂ < 0.

Test statistic

$$Z = \frac{P_1 - P_2}{\sqrt{P(1-P)\left (\frac{1}{n_1} + \frac{1}{n_2}\right)}}$$

where P is the proportion in the merged sample of size n₁ + n₂. Z ∈ N(0,1) and p-value can be computed using the standard normal distribution.

Two samples, mean

A two sample test of means is used to determine if two population means are equal.

Two independent samples are collected (one from each population) and the means are compared. Can for example be used to determine if a treatment group is different compared to a control group, in terms of the mean of a property of interest.

The null hypothesis;

H₀ : μ₂ = μ₁ The alternative hypothesis can either be two sided

H₁ : μ₂ ≠ μ₁ or one sided

H₁ : μ₂ > μ₁ or

H₁ : μ₂ < μ₁

Assume that observations from both populations are normally distributed;

$$ \begin{aligned} X_1 \sim N(\mu_1, \sigma_1) \\ X_2 \sim N(\mu_2, \sigma_2) \end{aligned} $$ Then it follows that the sample means will also be normally distributed;

$$ \begin{aligned} \bar X_1 \sim N(\mu_1, \sigma_1/\sqrt{n_1}) \\ \bar X_2 \sim N(\mu_2, \sigma_2/\sqrt{n_2}) \end{aligned} $$

The mean difference D = X̄₂ − X̄₁ is thus also normally distributed:

$$D = \bar X_2 - \bar X_1 = N\left(\mu_2-\mu_1, \sqrt{\frac{\sigma_2^2}{n_2} + \frac{\sigma_1^2}{n_1}}\right)$$

If H₀ is true: $$D = \bar X_2 - \bar X_1 = N\left(0, \sqrt{\frac{\sigma_2^2}{n_2} + \frac{\sigma_1^2}{n_1}}\right)$$

The test statistic: $$Z = \frac{\bar X_2 - \bar X_1}{\sqrt{\frac{\sigma_2^2}{n_2} + \frac{\sigma_1^2}{n_1}}}$$ is standard normal, i.e. Z ∼ N(0,1).

However, note that the test statistic require the standard deviations σ₁ and σ₂ to be known.

What if the population standard deviations are not known?

If the sample sizes are large, we can replace the known standard deviations with our sample standard deviations and according to the central limit theorem assume that

$$Z = \frac{\bar X_2 - \bar X_1}{\sqrt{\frac{s_2^2}{n_2} + \frac{s_1^2}{n_1}}} \sim N(0,1)$$

and proceed as before.

For small sample sizes the test statistic will be t-distributed.

$$t = \frac{\bar X_2 - \bar X_1}{\sqrt{\frac{s_2^2}{n_2} + \frac{s_1^2}{n_1}}}$$

For small sample sizes we can use Student’s t-test, which requires us to assume that X₁ and X₂ both are normally distributed and have equal variances. With these assumptions we can compute the pooled variance

$$ s_p^2 = \frac{(n_1-1)s_1^2 + (n_2-1)s_2^2}{n_1+n_2-2} $$

and the test statistic

$$t = \frac{\bar X_1 - \bar X_2}{\sqrt{s_p^2(\frac{1}{n_1} + \frac{1}{n_2})}}$$

t is t-distributed with n₁ + n₂ − 2 degrees of freedom.

The t-test is implemented in R, e.g. in the function t.test in the R-package stats, both Student’s t-test with equal variances and Welch’s t-test with unequal variances.

Variance

The test of equal variance in two groups is based on the null hypothesis

H₀ : σ₁² = σ₂²

If the two samples both come from two populations with normal distributions, the sample variances

$$S_1^2 = \frac{1}{n_1-1} \sum_{i=1}^{n_1} (X_{1i}-\bar X_1)^2\\ S_2^2 = \frac{1}{n_2-1} \sum_{i=1}^{n_2} (X_{2i}-\bar X_2)^2$$

It can be shown that $\frac{(n_1-1)S_1^2}{\sigma_1^2} \sim \chi^2(n_1-1)$ and $\frac{(n_2-1)S_2^2}{\sigma_2^2} \sim \chi^2(n_2-1)$.

Hence, the test statistic for comparing the variances of two groups

$$F = \frac{S_1^2}{S_2^2}$$ is F-distributed with n₁ − 1 and n₂ − 1 degrees of freedom.

In R a test of equal variances can be performed using the function var.test.

Exercises: Hypothesis tests, parametric

Does high fat diet increase body weight in mice?

1. Assume equal variances.
2. Don’t assume equal variances.

Solutions: Hypothesis tests, parametric

@ref(exr:Hb)

x <- c(154, 140, 147, 162, 172)
## sample mean
m <- mean(x)
## standard error of mean
SE <- sd(x)/sqrt(5)
## Observed value of test statistic
tobs <- (m-140)/SE
## P(t>tobs)
pt(tobs, df=4, lower.tail=FALSE)

## [1] 0.0277

## p = P(t>tobs) + P(t<-tobs) = 2 * P(t>tobs)
p <- 2*pt(tobs, df=4, lower.tail=FALSE)
p

## [1] 0.05541

@ref(exr:mice12) a)

# Student's t-test with pooled variances
t.test(xHF, xN, var.equal=TRUE, alternative="greater")

## 
##  Two Sample t-test
## 
## data:  xHF and xN
## t = 1, df = 22, p-value = 0.2
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  -1.469    Inf
## sample estimates:
## mean of x mean of y 
##     28.00     25.83

# Unequal variances with Welch approximation to the degrees of freedom (the default)
t.test(xHF, xN, var.equal=FALSE, alternative="greater")

## 
##  Welch Two Sample t-test
## 
## data:  xHF and xN
## t = 1, df = 20, p-value = 0.2
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  -1.486    Inf
## sample estimates:
## mean of x mean of y 
##     28.00     25.83

Point and interval estimates

As seen in previous chapter, the sample proportion or mean is an unbiased estimate of the population values. When we only have a sample, the sample estimate will be our best guess of the population value, but it will not be without error.

If we are interested in how large proportion of the Uppsala population is allergic to pollen, we can investigate this by studying a random sample. Randomly select 100 persons in Uppsala. It is important to actually sample randomly, ideally every individual should have the same probability of being sampled.

In our sample, we observe that 42 of the 100 has a pollen allergy. Hence, the observed sample proportion is p = 0.42.

Based on this observation our point estimate of the Uppsla popultation proportion π is π ≈ p = 0.42. We know that there is a certain uncertainty in this measurement, if the experiment is repeated we would select 100 other persons and our point estimate would be slightly different.

Bootstrap interval

Using bootstrap we can sample with replacement from our sample to estimate the uncertainty.

Bootstrap is to use the data we have (our sample) and sample repeatedly with replacement from this data.

Put the entire sample in an urn!

(#fig:pollenurn42)An urn model with 42 allergy (black) and 58 non-allergy (white). The black balls represent allergic and the white balls non-allergic.

Sample from the urn with replacement to compute the bootstrap distribution.

Using the bootstrap distribution the uncertainty of our estimate of π can be estimated.

The 95% bootstrap interval is [0.32, 0.52].

The bootstrap is very useful if you do not know the distribution of our sampled propery. But in our example we actually do.

Confidence interval

A confidence interval is a type of interval estimate associated with a confidence level.

Confidence interval of proportions

Remember that we can use the central limit theorem to show that

$$P \sim N\left(\pi, SE\right) \iff P \sim \left(\pi, \sqrt{\frac{\pi(1-\pi)}{n}}\right)$$

It follows that

$$Z = \frac{P - \pi}{SE} \sim N(0,1)$$ Based on what we know of the standard normal distribution, we can compute an interval around the population property π such that the probability that a sample property p fall within this interval is 1 − α.

$$P(-z < \frac{P - \pi}{SE} < z) = 1 - \alpha$$ For a 95% confidence interval z=1.96 (from a table of the standard normal distribution). Other confidence levels of interest include 90% (z=1.64) and 99% (z=2.58).

$$P\left(-z < \frac{P-\pi}{SE}<z\right) = \left(-z < Z <z\right) = 1-\alpha$$ We can rewrite this to

P(π−zSE<P<π+zSE) = 1 − α in other words sample fraction p will fall between π ± 1.96SE with 95% probability.

The equation can also be rewritten to P(P−zSE<π<P+zSE) = 1 − α The observed confidence interval is what we get when we replace the random variable P with our observed fraction,

p − zSE < π < p + zSE $$\pi = p \pm z SE = p \pm z \sqrt{\frac{p(1-p)}{n}}$$ The 95% confidence interval $$\pi = p \pm 1.96 \sqrt{\frac{p(1-p)}{n}}$$

A 95% confidence interval will have 95% chance to cover the true value.

Back to our example of proportion pollen allergic in Uppsala. p = 0.42 and $SE=\sqrt{\frac{p(1-p)}{n}} = 0.05$.

Hence, the 95% confidence interval is π = 0.42 ± 1.96 * 0.05 = 0.42 ± 0.092 or (0.42−0.092,0.42+0.092) = (0.32,0.52)

Confidence interval of mean

The confidence interval of mean can be derived similarly.

The mean of a sample of n independent and identically normal distributed observations X_i is normally distributed;

$$\bar X \sim N(\mu, \frac{\sigma}{\sqrt{n}})$$

If σ is unknown the statistic

$$\frac{\bar X - \mu}{\frac{\sigma}{\sqrt{n}}} \sim t(n-1)$$ is t-distributed with n − 1 degrees of freedom.

It follows that

$$ \begin{aligned} P\left(-t < \frac{\bar X - \mu}{\frac{\sigma}{\sqrt{n}}} < t\right) = 1 - \alpha \iff \\ P\left(\bar X - t \frac{\sigma}{\sqrt{n}} < \mu < \bar X + t \frac{\sigma}{\sqrt{n}}\right) = 1 - \alpha \end{aligned} $$

The confidence interval with confidence level 1 − α is thus;

$$\mu = \bar x \pm t \frac{s}{\sqrt{n}}$$

For a 95% confidence interval and n = 5 t is 2.78.

The t values for different values of α and degrees of freedom are tabulated and can be computed in R using the function qt.

n=5
alpha = 0.05
## t value
qt(1-alpha/2, df=n-1)

## [1] 2.776

Exercises: Interval estimates

Solutions

@ref(exr:CIHb)

obs <- c(145, 165, 134, 167, 158, 176, 156, 189, 143, 123)
mboot <- replicate(10000, {
  x <- sample(obs, size=10, replace=TRUE)
  mean(x)
})
hist(mboot)

## 05% confidence interval
quantile(mboot, c(0.025, 0.975))

##  2.5% 97.5% 
## 143.9 167.2

(m <- mean(obs))

## [1] 155.6

(v <- var(obs))

## [1] 395.2

(s <- sd(obs))

## [1] 19.88

5. Assume that Hb is normally distributed and compute a 95% confidence interval for the mean Hb value.

n <- length(obs)
t <- qt(0.975, df=9)
##95% confidence interval 
c(m - t*s/sqrt(n), m + t*s/sqrt(n))

## [1] 141.4 169.8

@ref(exr:CIprop)

1. Calculate a 90% confidence interval instead. Or sample more people than 100.

2. Change the z number,

p ± z * SE

For a 90% confidence interval use z=1.64

p <- 0.42
n <- 100
SE <- sqrt(p*(1-p)/n)
z <- qnorm(0.95)
c(p - z*SE, p + z*SE)

## [1] 0.3388 0.5012

3. or a 99% confidence interval?

z <- qnorm(0.995)
c(p - z*SE, p + z*SE)

## [1] 0.2929 0.5471

@(ref:CIscale)

The measured weight is a random variable X ∼ N(μ,σ). You know that σ = 2.3, μ is the weight of the smaple.

1. Compute a 95% confidence interval of the sample weight

## 95% ciónfidence interval
m <- 42
sigma <- 2.3
n <- 10
z <- qnorm(0.975)
c(m - z*sigma/sqrt(10), m + z*sigma/sqrt(10))

## [1] 40.57 43.43

2. Compute a 90% confidence interval of the sample weight

z <- qnorm(0.95)
c(m - z*sigma/sqrt(10), m + z*sigma/sqrt(10))

## [1] 40.8 43.2

@(ref:BMCsmokers)

You observe 150 students at BMC of which 25 are smokers. Compute a 95% confidence interval for the proportion of smokers among BMC students.

p <- 25/150
n <- 150
z <-qnorm(0.975)
SE <- sqrt(p*(1-p)/n)
## 95% CI
c(p - z*SE, p + z*SE)

## [1] 0.1070 0.2263

Multiple testing

Error types

Significance level

P(reject  H₀|H₀ is true) = P(type I error) = α

Statistical power

P(reject H₀|H₁ is true) = P(reject H₀|H₀ is false) = 1 − P(type II error)

Perform one test:

• P(One type I error) = α
• P(No type I error) = 1 − α

Perform m independent tests:

• P(No type I errors in m tests) = (1−α)^m
• P(At least one type I error in m tests) = 1 − (1−α)^m

• FWER: family-wise error rate, probability of one or more false positive, e.g. Bonferroni, Holm
• FDR: false discovery rate, proportion of false positives among hits, e.g. Benjamini-Hochberg, Storey

Bonferroni correction

To achieve a family-wise error rate of  ≤ α when performing m tests, declare significance and reject the null hypothesis for any test with p ≤ α/m.

Objections: too conservative

Benjamini-Hochbergs FDR

The false discovery rate is the proportion of false positives among ‘hits’, i.e. $\frac{FP}{TP+FP}$.

Benjamini-Hochberg’s method control the FDR level, γ, when performing m independent tests, as follows:

1. Sort the p-values p₁ ≤ p₂ ≤ … ≤ p_m.
2. Find the maximum j such that $p_j \leq \gamma \frac{j}{m}$.
3. Declare significance for all tests 1, 2, …, j.

‘Adjusted’ p-values

Sometimes an adjusted significance threshold is not reported, but instead ‘adjusted’ p-values are reported.

• Using Bonferroni’s method the ‘adjusted’ p-values are:

  p̃_i = min (mp_i,1).

A feature’s adjusted p-value represents the smallest FWER at which the null hypothesis will be rejected, i.e. the feature will be deemed significant.

• Benjamini-Hochberg’s ‘adjusted’ p-values are called q-values:

  $q_i = \min(\frac{m}{i} p_i, 1)$

  A feature’s q-value can be interpreted as the lowest FDR at which the corresponding null hypothesis will be rejected, i.e. the feature will be deemed significant.

Example, 10000 independent tests (e.g. genes)